A synthesis of dryland restoration techniques.

Purpose

To quantitatively examine the efficacy of vegetation restoration in drylands globally.

Questions

1. What is the global extent of research that directly examined restoration of drylands?
   
2. What were the common measures?
   
3. Is the restoration of vegetation a common and primary focus?
   
4. How frequently does the restoration measure outcomes beyond the focal species?
   
5. What were the primary restoration goals as reported by primary authors?
   
6. How much variation was there in the techniques tested and how long were experiments monitored and tested?
   
7. How relatively effective were the techniques?

Step 1. Search

Study-level viz to document patterns in exclusions primarily and the relatie frequenices, at the study level, of major categories of evidence.

#study data####
library(tidyverse)
studies <- read_csv("data/studies.csv")
studies

## # A tibble: 278 x 18
##       ID title technique data  region exclude rationale observations
##    <dbl> <chr> <chr>     <chr> <chr>  <chr>   <chr>     <chr>       
##  1   152 Shor~ seeding,~ expe~ Africa no      <NA>      <NA>        
##  2   180 Rest~ chemical~ App.~ Africa no      <NA>      <NA>        
##  3   229 Infl~ soil see~ fiel~ Africa no      <NA>      <NA>        
##  4   230 Acti~ planting  fiel~ Africa no      <NA>      <NA>        
##  5   255 The ~ grazing ~ fiel~ Africa no      <NA>      <NA>        
##  6   262 Reve~ seeding,~ eper~ Africa no      <NA>      <NA>        
##  7   263 The ~ phytogen~ fiel~ Africa no      <NA>      <NA>        
##  8   264 Eval~ seeding,~ fiel~ Africa no      <NA>      <NA>        
##  9   271 Patc~ natural ~ fiel~ Africa no      <NA>      <NA>        
## 10     4 Fact~ natural ~ App.~ Africa no      <NA>      <NA>        
## # ... with 268 more rows, and 10 more variables: disturbance <chr>,
## #   system <chr>, goal <chr>, intervention <chr>, paradigm <chr>,
## #   grazing <chr>, hypothesis <chr>, soil <chr>, benchmark <chr>,
## #   notes <chr>

#quick look at rationale needed
exclusions <- studies %>%
  filter(exclude == "yes")

#quick look at studies with paradigms
evidence <- studies %>%
  filter(exclude == "no")

#library(skimr)
#skim(evidence)

#study-level viz#####
#exclusions
#ggplot(exclusions, aes(rationale, fill = region)) +
 # geom_bar() +
 # coord_flip() +
 # labs(x = "rational for exclusion", y = "frequency") +
 # scale_fill_brewer(palette = "Paired")

ggplot(evidence, aes(disturbance, fill = paradigm)) +
  geom_bar(na.rm = TRUE) +
  coord_flip() +
  scale_fill_brewer(palette = "Paired") +
  labs(y = "frequency")

#ggplot(evidence, aes(region, fill = paradigm)) +
 # geom_bar(na.rm = TRUE) +
 # coord_flip() +
 # scale_fill_brewer(palette = "Paired") +
 # labs(y = "frequency")

#ggplot(evidence, aes(data, fill = paradigm)) +
 # geom_bar(na.rm = TRUE) +
  #coord_flip() +
  #scale_fill_brewer(palette = "Paired") +
  #labs(y = "frequency")

#ggplot(evidence, aes(system, fill = paradigm)) +
 # geom_bar(na.rm = TRUE) +
 # coord_flip() +
 # scale_fill_brewer(palette = "Paired") +
 # labs(y = "frequency")

#ggplot(evidence, aes(goal, fill = paradigm)) +
  #geom_bar(na.rm = TRUE) +
  #coord_flip() +
  #scale_fill_brewer(palette = "Paired") +
  #labs(x = "outcome", y = "frequency")

#step 1 models####
#paradigm
derived.evidence <- evidence %>%
  group_by(technique, data, region, disturbance, goal, paradigm) %>% summarise(n = n())

#active-passive split
m <- glm(n~paradigm, family = poisson, derived.evidence)
anova(m, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: n
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                       167     9.9147         
## paradigm  1 0.045115       166     9.8696   0.8318

#region
m1 <- glm(n~paradigm*region, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: n
## 
## Terms added sequentially (first to last)
## 
## 
##                 Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                              167     9.9147         
## paradigm         1 0.045115       166     9.8696   0.8318
## region           6 0.301367       160     9.5682   0.9995
## paradigm:region  6 0.213627       154     9.3546   0.9998

#outcome
m2 <- glm(n~paradigm*goal, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m2, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: n
## 
## Terms added sequentially (first to last)
## 
## 
##               Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                            167     9.9147         
## paradigm       1 0.045115       166     9.8696   0.8318
## goal           6 0.240941       160     9.6287   0.9997
## paradigm:goal  4 0.301480       156     9.3272   0.9897

#even split between active and passive evidence by all key categories

Step 2. Sort

A summary of sort process using PRISMA.

library(PRISMAstatement)
prisma(found = 1504,
       found_other = 5,
       no_dupes = 1039, 
       screened = 1039, 
       screen_exclusions = 861, 
       full_text = 178,
       full_text_exclusions = 101, 
       qualitative = 77, 
       quantitative = 66,
       width = 800, height = 800)

Step 3. Synthesize

Check data and calculate necessary measures.

#all data
#data <- read_csv("data/data.csv")
#data <- data %>%
  #mutate(lrr = log(mean.t/mean.c), rii = ((mean.t-mean.c)/(mean.t + mean.c)), var.es = ((sd.t^2/n.t*mean.t^2) + (sd.c^2/n.c*mean.c^2)))

#other effect size estimates
#library(compute.es)
#data <- data %>%
  #mutate(d=mes(mean.t, mean.c, sd.t, sd.c, n.t, n.c, level = 95, #cer = 0.2, dig = 2, , id = ID, data = data))
#check metafor

#data from ag and grazing studies that examined restoration in drylands
data <- read_csv("data/data.csv")
mydata<- data %>% 
 filter(disturbance %in% c("agriculture","grazing")) %>% 
filter(!notes %in% "couldnt extract data")

#write.csv(mydata, file="mydata.csv")

mydata <- mydata %>%
  mutate(lrr = log(mean.t/mean.c), rii = ((mean.t-mean.c)/(mean.t + mean.c)), var.es = ((sd.t^2/(n.t*mean.t^2)) + (sd.c^2/(n.c*mean.c^2)))) #use only lrr now

#mydata<- mydata %>% 
 # mutate(aridity.range = cut(aridity.index, breaks=c(0,5,10,20,25,69.92), labels=c("hyperarid","arid", #"semiarid", "moderately arid", "slightly humid"))) #categorization of aridity.index values, according to de Martone

Step 4. Summarize

Explore summary level data of all data. Explore aggregation levels that support the most reasonable data structure and minimize non-independence issues.

#evidence map####
require(maps)
world<-map_data("world")
map<-ggplot() + geom_polygon(data=world, fill="gray50", aes(x=long, y=lat, group=group))
#map + geom_point(data=paperdata, aes(x=long, y=lat)) +
  #labs(x = "longitude", y = "latitude") #render a literal map, i.e. evidence map, of where we study the niche in deserts globally

#add in levels and color code points on map####
#map of 178 articles
map + geom_point(data=data, aes(x=long, y=lat, color = paradigm)) + 
  scale_color_brewer(palette = "Paired") +
  labs(x = "longitude", y = "latitude", color = "")

#map of selected articles, agriculture and grazing disturbances
map + geom_point(data=mydata, aes(x=long, y=lat, color = paradigm)) + 
  scale_color_brewer(palette = "Paired") +
  labs(x = "longitude", y = "latitude", color = "")

#aggregation####
se <- function(x){
  sd(x)/sqrt(length(x))
}

data.simple <- mydata %>%
  group_by(study.ID, paradigm, technique, measure.success) %>%
  summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))

main.data <- mydata %>%
  group_by(study.ID, paradigm, intervention, outcome) %>%
  summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))


#EDA data####
simple.data <- mydata %>% group_by(study.ID, paradigm, technique, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
simple.data <- na.omit(simple.data)

parad.data <- mydata %>% group_by(study.ID, paradigm) %>% summarise(mean.rii = mean(rii), error = se(rii))
parad.data <- na.omit(parad.data)

tech.data <- mydata %>% group_by(study.ID, technique) %>% summarise(mean.rii = mean(rii), error = se(rii))
tech.data <- na.omit(tech.data)

success.data <- mydata %>% group_by(study.ID, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
success.data <- na.omit(success.data)


#active
active <- mydata %>%
  filter(paradigm == "active")

#viz for aggregation####
disturbance.data <- mydata %>% 
  group_by(study.ID,disturbance, paradigm) %>% 
  count()

disturbance.data2 <- disturbance.data %>% 
  group_by(disturbance,paradigm) %>% 
  count()

ggplot(na.omit(disturbance.data2), aes(disturbance,nn, fill=paradigm))+ 
  geom_bar(stat = "identity") + 
  coord_flip(ylim=0:44) + 
  scale_fill_brewer(palette = "Blues")

intervention.data <- active %>% 
  group_by(study.ID,intervention, paradigm) %>% 
  count()

intervention.data2 <- intervention.data %>% 
  group_by(intervention,paradigm) %>% 
  count()

#ggplot(na.omit(intervention.data2), aes(intervention,nn, fill=paradigm)) +
  #geom_bar(stat = "identity") + 
  #coord_flip(ylim=0:44) + 
  #scale_fill_brewer(palette = "Blues")

technique.data <- mydata %>% 
  group_by(study.ID,technique, paradigm) %>% 
  count()

technique.data2 <- technique.data %>% 
  group_by(technique,paradigm) %>% 
  count()

technique.data3 <- technique.data2 %>% 
  group_by(technique) %>% 
  count()

aridity<- mydata %>% 
  group_by(continent) %>% summarize(mean.aridity=mean(aridity.index))
aridity  

## # A tibble: 6 x 2
##   continent     mean.aridity
##   <chr>                <dbl>
## 1 Africa                9.04
## 2 Asia                 20.7 
## 3 Europe               26.8 
## 4 North America        19.6 
## 5 Oceania              16.4 
## 6 South America        22.9

#table(mydata$aridity.index)


ggplot(na.omit(data.simple), aes(technique, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired")

ggplot(na.omit(data.simple), aes(measure.success, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired")

#better
ggplot(main.data, aes(intervention, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired") +
  labs(fill = "")

ggplot(main.data, aes(outcome, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired") +
  labs(fill = "")

Step 5a. EDA

Exploratory data analyses to understand data and QA/QC using Rii.

Step 5b. Models

Meta and conventional statistical models to explore relative efficacy.

Step 5. Synthesis stats

#p-value meta
#library(metap)
#all data for metas but cleaned for na's
#all data for meta cleaned
mdata.all <- mydata %>%
  filter(!is.na(lrr)) %>%
  filter(!is.na(var.es)) %>%
  filter(!is.na(n.t)) %>%
  filter(!is.na(p)) %>%
  filter(!is.na(intervention)) %>%
  filter(is.finite(lrr)) %>%
  filter(!is.na(exp.length)) %>%
  filter(!is.na(MAP)) %>%
  filter(!is.na(aridity.index)) #%>% 
  #filter(!is.infinite(var.es))


summary(mdata.all$var.es)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.00000  0.00000  0.00076  0.15430  0.01441 39.06266

plot(density(mdata.all$var.es))

top_n(mdata.all, -10, var.es)

## # A tibble: 10 x 50
##    study.ID    ID publication.year data.year exp.year exp.length
##       <dbl> <dbl>            <dbl>     <dbl>    <dbl>      <dbl>
##  1      111  111.             2016      2003     2003        132
##  2      111  111.             2016      2003     2003        132
##  3      111  111.             2016      2003     2003        132
##  4      111  111.             2016      2003     2003        132
##  5      111  111.             2016      2003     2003        132
##  6      111  111.             2016      2003     2003        132
##  7      111  111.             2016      2003     2003        132
##  8      111  111.             2016      2003     2003        132
##  9      111  111.             2016      2003     2003        132
## 10      111  111.             2016      2003     2003        132
## # ... with 44 more variables: disturbance <chr>, focus <chr>,
## #   technique <chr>, intervention <chr>, paradigm <chr>, hypothesis <chr>,
## #   pathway <chr>, plant.species <chr>, target.plant <chr>,
## #   measure.success <chr>, outcome <chr>, measured.factor <chr>,
## #   treatment <chr>, control <chr>, unit <chr>, Nsites <dbl>, n.t <dbl>,
## #   n.c <dbl>, ntotalsamples <dbl>, mean.t <dbl>, mean.c <dbl>,
## #   sd.t <dbl>, sd.c <dbl>, se.t <dbl>, se.c <dbl>, p <dbl>, df <dbl>,
## #   measure.dispersion <chr>, lat <dbl>, long <dbl>, continent <chr>,
## #   country <chr>, ecosystem <chr>, elevation <dbl>, MAP <dbl>, MAT <dbl>,
## #   aridity.index <dbl>, `potential.evaporation (mm)` <dbl>,
## #   grazing <chr>, soil <chr>, notes <chr>, lrr <dbl>, rii <dbl>,
## #   var.es <dbl>

top_n(mdata.all, 5, var.es)

## # A tibble: 5 x 50
##   study.ID    ID publication.year data.year exp.year exp.length disturbance
##      <dbl> <dbl>            <dbl>     <dbl>    <dbl>      <dbl> <chr>      
## 1       11  11.7             2015      2007     1956        600 agriculture
## 2       11  11.7             2015      2007     1956        600 agriculture
## 3       11  11.7             2015      2007     1956        600 agriculture
## 4       11  11.7             2015      2007     1956        600 agriculture
## 5      147 147.              2012      2005     2005         36 agriculture
## # ... with 43 more variables: focus <chr>, technique <chr>,
## #   intervention <chr>, paradigm <chr>, hypothesis <chr>, pathway <chr>,
## #   plant.species <chr>, target.plant <chr>, measure.success <chr>,
## #   outcome <chr>, measured.factor <chr>, treatment <chr>, control <chr>,
## #   unit <chr>, Nsites <dbl>, n.t <dbl>, n.c <dbl>, ntotalsamples <dbl>,
## #   mean.t <dbl>, mean.c <dbl>, sd.t <dbl>, sd.c <dbl>, se.t <dbl>,
## #   se.c <dbl>, p <dbl>, df <dbl>, measure.dispersion <chr>, lat <dbl>,
## #   long <dbl>, continent <chr>, country <chr>, ecosystem <chr>,
## #   elevation <dbl>, MAP <dbl>, MAT <dbl>, aridity.index <dbl>,
## #   `potential.evaporation (mm)` <dbl>, grazing <chr>, soil <chr>,
## #   notes <chr>, lrr <dbl>, rii <dbl>, var.es <dbl>

library(metafor)

#metadata<-escalc(measure="ROM",m1i=mean.t,m2i=mean.c,sd1i=sd.t,sd2i=sd.c,n1i=n.t,n2i=n.c, #data=mydata,var.names=c("LRR","LRR_var"),digits=4)

#metadata <- metadata %>%
#  filter(!is.na(LRR)) %>%
#  filter(!is.na(LRR_var)) %>%
#  filter(!is.na(n.t)) %>%
#  filter(!is.na(p)) %>%
#  filter(!is.na(intervention)) %>%
#  filter(is.finite(lrr)) %>%
#  filter(!is.na(exp.length)) %>% 
#  filter(!is.na(aridity.index))

#summary(metadata$LRR_var)
#plot(density(metadata$var.es))
#top_n(metadata, -10, LRR_var)
#top_n(metadata, 5, LRR_var)

max(mdata.all$var.es)/min(mdata.all$var.es)

## [1] 4.800933e+14

#active only data for meta cleaned up
mdata <- mydata %>%
  filter(paradigm == "active") %>%
  filter(!is.na(lrr)) %>%
  filter(!is.na(var.es)) %>%
  filter(!is.na(n.t)) %>%
  filter(!is.na(p)) %>%
  filter(!is.na(intervention)) %>%
  filter(is.finite(lrr)) %>%
  filter(!is.na(exp.length))

#aggregated data for metas var estimated with central tendency
#note - could also bootstrap mean variance here instead of arithematic mean
simple.mdata <- mdata %>%
  group_by(intervention) %>%
  summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))

simple.mdata.2 <- mdata %>%
  group_by(intervention, outcome) %>%
  summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))  

simple.mdata.var <- mdata %>%
  group_by(intervention) %>%
  summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))

simple.mdata2.var <- mdata %>%
  group_by(intervention, outcome) %>%
  summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))

#metas with p-values####
#schweder(mdata$p)
#sumz(p, data = mdata)
#mdata %>%
  #split(.$intervention) %>%
  #purrr::map(~sumz(.$p)) 
#sumlog(mdata$p)

#metas with meta package on effect size measures####
library(meta)
#all data non-aggregated
#m <- metagen(lrr, var.es, studlab = ID, byvar = intervention, data = mdata)
#summary(m)
#funnel(m)
#metabias(m)
#forest(m, sortvar = intervention)

#t-tests if different from 0

tmu <- function(x){t.test(x, mu = 0, paired = FALSE, var.equal=FALSE, conf.level = 0.95)}

mdata.all %>%
  split(.$paradigm) %>%
  purrr::map(~tmu(.$lrr)) #note this uses arithmetic means not estimated means from random effect models

## $active
## 
##  One Sample t-test
## 
## data:  x
## t = 7.6083, df = 1101, p-value = 5.943e-14
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.2069479 0.3507825
## sample estimates:
## mean of x 
## 0.2788652 
## 
## 
## $passive
## 
##  One Sample t-test
## 
## data:  x
## t = -7.5438, df = 357, p-value = 3.824e-13
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.4408271 -0.2585133
## sample estimates:
##  mean of x 
## -0.3496702

#Propose we only use random effect models because of the diversity of studies, differences in the methods and samples that may introduce heterogeneity
m1 <- metagen(lrr, var.es, studlab = ID, byvar = paradigm, comb.fixed=FALSE, data = mdata.all)
summary(m1) #active is positive and passive is negative so should not mix

## Number of studies combined: k = 1460
## 
##                                       95%-CI     z  p-value
## Random effects model 0.0766 [0.0654; 0.0879] 13.38 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 184527270028.74 [184527270027.80; 184527270029.68]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 184583923516.49 [184583923515.55; 184583923517.44]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                              Q d.f. p-value
##  49679407227636228416848646.00 1459       0
## 
## Results for subgroups (random effects model):
##                       k                     95%-CI
## paradigm = active  1102  0.2184 [ 0.2055;  0.2314]
## paradigm = passive  358 -0.3413 [-0.3753; -0.3073]
##                                                Q  tau^2    I^2
## paradigm = active  49671693556322550580400868.00 0.0450 100.0%
## paradigm = passive     4152232320954931871684.00 0.1047 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   911.23    1 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#metabias(m)
#should do a t-test right now of paradigm is diff from 0
metareg(m1,aridity.index)

## 
## Mixed-Effects Model (k = 1460; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value):             0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   33269628721436362276444.00
## R^2 (amount of heterogeneity accounted for):            2.20%
## 
## Test for Residual Heterogeneity:
## QE(df = 1458) = 48507118675854218692260048.0000, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 2845.3011, p-val < .0001
## 
## Model Results:
## 
##                estimate      se      zval    pval    ci.lb    ci.ub 
## intrcpt          0.7855  0.0144   54.3680  <.0001   0.7571   0.8138  *** 
## aridity.index   -0.0354  0.0007  -53.3414  <.0001  -0.0367  -0.0341  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#interventions
m2 <- metagen(lrr, var.es, studlab = ID, byvar = intervention, comb.fixed=FALSE, subset = paradigm == "active", data = mdata.all)
summary(m2)

## Number of studies combined: k = 1102
## 
##                                       95%-CI     z  p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 210054235083.87 [210054235082.78; 210054235084.96]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                              Q d.f. p-value
##  49671693556322550580400868.00 1101       0
## 
## Results for subgroups (random effects model):
##                                 k                  95%-CI
## intervention = vegetation     779 0.1845 [0.1694; 0.1996]
## intervention = soil           248 0.3128 [0.2990; 0.3265]
## intervention = water addition  75 0.6409 [0.5539; 0.7279]
##                                                           Q  tau^2    I^2
## intervention = vegetation     48393423300974332080088226.00 0.0443 100.0%
## intervention = soil              97513761686148203152860.00 0.0111 100.0%
## intervention = water addition                      31132.18 0.1047  99.8%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   226.58    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#metabias(m)

# Meta regression for aridity index (de Martonne 1927) and length of experiments (in months)
metareg(m2, ~ aridity.index)

## 
## Mixed-Effects Model (k = 1102; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value):             0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   44082670438161540188602.00
## R^2 (amount of heterogeneity accounted for):            2.22%
## 
## Test for Residual Heterogeneity:
## QE(df = 1100) = 48490937481977699572460620.0000, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 1474.0362, p-val < .0001
## 
## Model Results:
## 
##                estimate      se      zval    pval    ci.lb    ci.ub 
## intrcpt          0.7703  0.0158   48.7709  <.0001   0.7393   0.8012  *** 
## aridity.index   -0.0289  0.0008  -38.3932  <.0001  -0.0304  -0.0274  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

metareg(m2, ~ exp.length )

## 
## Mixed-Effects Model (k = 1102; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value):             0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   44082669939996730852248.00
## R^2 (amount of heterogeneity accounted for):            2.22%
## 
## Test for Residual Heterogeneity:
## QE(df = 1100) = 48490936933996402964282222.0000, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 961.4837, p-val < .0001
## 
## Model Results:
## 
##             estimate      se     zval    pval   ci.lb   ci.ub 
## intrcpt       0.0499  0.0085   5.8663  <.0001  0.0332  0.0666  *** 
## exp.length    0.0014  0.0000  31.0078  <.0001  0.0013  0.0015  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m3 <- metagen(lrr, var.es, studlab = ID, byvar = intervention, subset = paradigm == "passive", comb.fixed=FALSE, data = mdata.all)
summary(m3)

## Number of studies combined: k = 358
## 
##                                          95%-CI      z  p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 3420000814.26 [3420000812.64; 3420000815.87]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                          Q d.f. p-value
##  4152232320954931871684.00  357       0
## 
## Results for subgroups (random effects model):
##                                    k                     95%-CI
## intervention = vegetation        125  0.2654 [ 0.2067;  0.3241]
## intervention = grazing exclusion  29  0.1351 [ 0.0270;  0.2431]
## intervention = soil              204 -0.7583 [-0.8196; -0.6970]
##                                                          Q  tau^2    I^2
## intervention = vegetation        4152209524073903423668.00 0.1047 100.0%
## intervention = grazing exclusion              238316232.18 0.0881 100.0%
## intervention = soil                   14453104123321616.00 0.1990 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   595.91    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#metabias(m)


#trying models by continents
m4 <- metagen(lrr, var.es, studlab = ID, byvar = continent, subset = paradigm == "active", comb.fixed=FALSE, data = mdata.all)
summary(m4)

## Number of studies combined: k = 1102
## 
##                                       95%-CI     z  p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 210341521299.85 [210341521298.76; 210341521300.95]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                              Q d.f. p-value
##  49671693556322550580400868.00 1101       0
## 
## Results for subgroups (random effects model):
##                             k                     95%-CI
## continent = Asia          652  0.2110 [ 0.1947;  0.2273]
## continent = Europe        253  0.3272 [ 0.3136;  0.3408]
## continent = North America 150  0.0531 [ 0.0245;  0.0816]
## continent = South America  19 -0.4895 [-0.6073; -0.3718]
## continent = Africa         16  1.6107 [ 1.3665;  1.8549]
## continent = Oceania        12 -0.0398 [-0.1834;  0.1037]
##                                                       Q  tau^2    I^2
## continent = Asia          48393423133721090038820246.00 0.0443 100.0%
## continent = Europe           97513761686148203152860.00 0.0111 100.0%
## continent = North America          23018426431807620.00 0.0209 100.0%
## continent = South America                      90607.11 0.0643 100.0%
## continent = Africa                   254005421361935.34 0.1885 100.0%
## continent = Oceania                             1485.24 0.0628  99.3%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   615.00    5 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

m5 <- metagen(lrr, var.es, studlab = ID, byvar = continent, subset = paradigm == "passive", comb.fixed=FALSE, data = mdata.all)
summary(m5)

## Number of studies combined: k = 358
## 
##                                          95%-CI      z  p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 3429673787.56 [3429673785.94; 3429673789.18]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                          Q d.f. p-value
##  4152232320954931871684.00  357       0
## 
## Results for subgroups (random effects model):
##                             k                     95%-CI
## continent = North America  24  0.1711 [ 0.1162;  0.2260]
## continent = Asia          284 -0.4924 [-0.5448; -0.4400]
## continent = Europe         21  0.0886 [-0.0387;  0.2159]
## continent = South America   4  0.7712 [-0.6724;  2.2148]
## continent = Africa         25  0.3817 [ 0.2547;  0.5087]
##                                                   Q  tau^2    I^2
## continent = North America                    709.76 0.0144  96.8%
## continent = Asia               14453104145902294.00 0.1990 100.0%
## continent = Europe                     238286159.01 0.0886 100.0%
## continent = South America      33314949487640844.00 2.1620 100.0%
## continent = Africa        4152172019980625248426.00 0.1047 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   373.50    4 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#use random effects mean and var estimate to plot

#do t-tests here too

#outcomes
m <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "active", comb.fixed=FALSE, data = mdata.all)
summary(m)

## Number of studies combined: k = 1102
## 
##                                       95%-CI     z  p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 210149866430.08 [210149866428.99; 210149866431.17]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                              Q d.f. p-value
##  49671693556322550580400868.00 1101       0
## 
## Results for subgroups (random effects model):
##                     k                     95%-CI
## outcome = soil    249  0.2204 [ 0.1558;  0.2849]
## outcome = plants  305  0.5071 [ 0.4936;  0.5206]
## outcome = animals  24 -0.1152 [-0.1155; -0.1148]
## outcome = habitat 524  0.0621 [ 0.0437;  0.0804]
##                                               Q   tau^2    I^2
## outcome = soil                35077220764051.67  0.2656 100.0%
## outcome = plants     97513760543782884872406.00  0.0111 100.0%
## outcome = animals                     541696.64 <0.0001 100.0%
## outcome = habitat 48393423303339003636244404.00  0.0443 100.0%
## 
## Test for subgroup differences (random effects model):
##                        Q d.f. p-value
## Between groups   8647.81    3       0
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#metabias(m)

m <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "passive", comb.fixed=FALSE, data = mdata.all)
summary(m)

## Number of studies combined: k = 358
## 
##                                          95%-CI      z  p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 3419999089.05 [3419999087.44; 3419999090.67]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                          Q d.f. p-value
##  4152232320954931871684.00  357       0
## 
## Results for subgroups (random effects model):
##                     k                     95%-CI                         Q
## outcome = habitat 104  0.1605 [ 0.0964;  0.2246] 4152172019980636258444.00
## outcome = plants   50  0.4438 [ 0.0345;  0.8532]      33314950066906688.00
## outcome = soil    204 -0.7583 [-0.8196; -0.6970]      14453104123321616.00
##                    tau^2    I^2
## outcome = habitat 0.1047 100.0%
## outcome = plants  2.1620 100.0%
## outcome = soil    0.1990 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   425.72    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#metabias(m)

#t-tests for outcomes diff from 0 but you can see using the 95% CI what is different

#check metafor for interactions?? in these big models or are we ok??

#brainstorm on how to explore?? techniques


#save but cut all this.
#m.study <- metagen(lrr, var.es, studlab = study.ID, byvar = intervention, data = mdata)
#summary(m.study)
#funnel(m)
#metabias(m)
#forest(m, sortvar = intervention)

#aggregated data
#m1 <- metagen(lrr, var.es, studlab = intervention, data = simple.mdata)
#summary(m1)
#funnel(m1)
#metabias(m1)
#forest(m, sortvar = intervention)

#different var estimate
#m2 <- metagen(lrr, var.es, studlab = intervention, data = simple.mdata.var)
#summary(m2)
#funnel(m2)
#metabias(m2)
#forest(m, sortvar = intervention)


#m3 <- metagen(lrr, var.es, studlab = intervention, byvar = outcome, data = simple.mdata.2)
#summary(m3)
#funnel(m)
#radial(m3)
#metabias(m2)
#forest(m, sortvar = intervention)

#m4 <- metagen(lrr, var.es, studlab = intervention, byvar = outcome, data = simple.mdata2.var)
#summary(m4)
#funnel(m)
#radial(m4)
#metabias(m2)
#forest(m, sortvar = intervention)

library(metafor)


# I´ve installed the 'devel' version of metafor (https://wviechtb.github.io/metafor/#installation)

#checking and comparing Lrr and sampling variances values obtained with escalc() function and by functions manually defined (lines 153-155)

#metadata<-escalc(measure="ROM",m1i=mean.t,m2i=mean.c,sd1i=sd.t,sd2i=sd.c,n1i=n.t,n2i=n.c, #data=mydata,var.names=c("LRR","LRR_var"),digits=4)


#metadata <- metadata %>%
#  filter(!is.na(LRR)) %>%
#  filter(!is.na(LRR_var)) %>%
#  filter(!is.na(n.t)) %>%
#  filter(!is.na(p)) %>%
#  filter(!is.na(intervention)) %>%
#  filter(is.finite(lrr)) %>%
#  filter(!is.na(exp.length)) %>% 
#  filter(!is.na(aridity.index))

#summary(metadata$LRR_var)
#plot(density(metadata$var.es))
#top_n(metadata, -10, LRR_var)
#top_n(metadata, 5, LRR_var)


#write.csv(metadata, "metadata.csv")

#mdata.veg <- mydata %>%
 # filter(paradigm == "active") %>%
  #filter(intervention=="vegetation") %>% #filtering by intervention
  #filter(!is.na(lrr)) %>%
  #filter(!is.na(var.es)) %>%
#  filter(!is.na(n.t)) %>%
#  filter(!is.na(p)) %>%
#  filter(!is.na(intervention)) %>%
#  filter(is.finite(lrr)) %>%
#  filter(!is.na(exp.length))

#mdata.mean<- mdata.all %>% 
#  group_by(study.ID, aridity.index, exp.length, paradigm, intervention, outcome) %>% 
#  summarize(mean_lrr= mean(lrr),
#  mean_var.es= mean(var.es)) #collapsing data to means

#mdata.mean<- mdata.mean %>%
#  filter(!is.infinite(mean_var.es)) #filter out outliers that are infinite


#res1 <- rma(lrr, var.es, mods= aridity.index, data=mdata.all, subset=intervention=="vegetation") 
#res2 <- rma(lrr, var.es, data=mdata.all, subset=intervention=="soil") error

mod.1<-rma(lrr,var.es, mods= ~ aridity.index, subset= paradigm == "active", digits=4,data=mdata.all)
mod.1

## 
## Mixed-Effects Model (k = 1102; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.8019 (SE = 0.0353)
## tau (square root of estimated tau^2 value):             0.8955
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   113092375770.11
## R^2 (amount of heterogeneity accounted for):            8.78%
## 
## Test for Residual Heterogeneity:
## QE(df = 1100) = 6499005541922.0801, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 112.0030, p-val < .0001
## 
## Model Results:
## 
##                estimate      se      zval    pval    ci.lb    ci.ub 
## intrcpt          0.8574  0.0651   13.1675  <.0001   0.7298   0.9851  *** 
## aridity.index   -0.0330  0.0031  -10.5831  <.0001  -0.0391  -0.0269  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mod.2<-rma.mv(yi=lrr, V=var.es, mods= aridity.index, method = "REML",  test="t", random = ~ 1 | ID, data=mdata.all, sparse=TRUE)

mod.3 <- rma(yi=lrr, vi=var.es, data = mdata.all)
summary(m3)

## Number of studies combined: k = 358
## 
##                                          95%-CI      z  p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 3420000814.26 [3420000812.64; 3420000815.87]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                          Q d.f. p-value
##  4152232320954931871684.00  357       0
## 
## Results for subgroups (random effects model):
##                                    k                     95%-CI
## intervention = vegetation        125  0.2654 [ 0.2067;  0.3241]
## intervention = grazing exclusion  29  0.1351 [ 0.0270;  0.2431]
## intervention = soil              204 -0.7583 [-0.8196; -0.6970]
##                                                          Q  tau^2    I^2
## intervention = vegetation        4152209524073903423668.00 0.1047 100.0%
## intervention = grazing exclusion              238316232.18 0.0881 100.0%
## intervention = soil                   14453104123321616.00 0.1990 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   595.91    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

mod.4 <- rma(lrr, var.es, mods= cbind (aridity.index, exp.length), data = mdata.all, subset = paradigm == "active")
summary(m4)

## Number of studies combined: k = 1102
## 
##                                       95%-CI     z  p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 210341521299.85 [210341521298.76; 210341521300.95]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                              Q d.f. p-value
##  49671693556322550580400868.00 1101       0
## 
## Results for subgroups (random effects model):
##                             k                     95%-CI
## continent = Asia          652  0.2110 [ 0.1947;  0.2273]
## continent = Europe        253  0.3272 [ 0.3136;  0.3408]
## continent = North America 150  0.0531 [ 0.0245;  0.0816]
## continent = South America  19 -0.4895 [-0.6073; -0.3718]
## continent = Africa         16  1.6107 [ 1.3665;  1.8549]
## continent = Oceania        12 -0.0398 [-0.1834;  0.1037]
##                                                       Q  tau^2    I^2
## continent = Asia          48393423133721090038820246.00 0.0443 100.0%
## continent = Europe           97513761686148203152860.00 0.0111 100.0%
## continent = North America          23018426431807620.00 0.0209 100.0%
## continent = South America                      90607.11 0.0643 100.0%
## continent = Africa                   254005421361935.34 0.1885 100.0%
## continent = Oceania                             1485.24 0.0628  99.3%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   615.00    5 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

mod.5 <- rma(lrr, var.es, mods=  ~intervention, data = mdata.all, subset = paradigm == "active")
summary(m5)

## Number of studies combined: k = 358
## 
##                                          95%-CI      z  p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Quantifying residual heterogeneity:
## H = 3429673787.56 [3429673785.94; 3429673789.18]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                          Q d.f. p-value
##  4152232320954931871684.00  357       0
## 
## Results for subgroups (random effects model):
##                             k                     95%-CI
## continent = North America  24  0.1711 [ 0.1162;  0.2260]
## continent = Asia          284 -0.4924 [-0.5448; -0.4400]
## continent = Europe         21  0.0886 [-0.0387;  0.2159]
## continent = South America   4  0.7712 [-0.6724;  2.2148]
## continent = Africa         25  0.3817 [ 0.2547;  0.5087]
##                                                   Q  tau^2    I^2
## continent = North America                    709.76 0.0144  96.8%
## continent = Asia               14453104145902294.00 0.1990 100.0%
## continent = Europe                     238286159.01 0.0886 100.0%
## continent = South America      33314949487640844.00 2.1620 100.0%
## continent = Africa        4152172019980625248426.00 0.1047 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   373.50    4 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

#forest(m1, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#forest(m2, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#forest(m3, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#forest(m4, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#t-tests for lrr diff from 0
#mdata %>%
  #split(.$intervention) %>%
  #purrr::map(~sumz(.$lrr)) 

#effect sizes plots
#need better ci estimates
#ggplot(simple.mdata, aes(intervention, lrr)) +
 # ylim(c(-2,2)) +
#  geom_point(position = position_dodge(width = 0.5)) + 
 # labs(x = "", y = "lrr") +
 # coord_flip() +
 # geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.05, position = position_dodge(width = 0.5)) +
 # geom_hline(yintercept = 0, colour="grey", linetype = "longdash")

#ggplot(simple.mdata.2, aes(intervention, lrr, color = outcome)) +
 # ylim(c(-2,2)) +
  #geom_point(position = position_dodge(width = 0.5)) + 
 # labs(x = "", y = "lrr", color = "") +
 # coord_flip() +
  #geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.05, position = #position_dodge(width = 0.5)) +
  #geom_hline(yintercept = 0, colour="grey", linetype = "longdash")

#ggplot(mdata, aes(lrr, color = intervention)) +
  #geom_freqpoly(binwidth = 0.5, size = 2) + 
  #xlim(c(-10, 10)) +
  #labs(x = "lrr", y = "frequency", color = "") +
  #geom_vline(xintercept = 0, colour="grey", linetype = "longdash") +
  #scale_color_brewer()

#ggplot(mdata, aes(lrr, fill = intervention)) +
  #geom_dotplot(binwidth = 1) + 
 # xlim(c(-10, 10)) +
 # labs(x = "lrr", y = "frequency", fill = "") +
 # geom_vline(xintercept = 0, colour="grey", linetype = "longdash") +
 # scale_fill_brewer()

Interpretations

1. Little to no replication of specific techniques.
   
2. Intervention and outcome classification useful and appropriate.
   
3. Active versus passive interesting and fascinating.
   
4. Use only random effect model statistics because of heterogeneity and conceptual focus of data and synthesis.
   
5. Active and passive are not the same - active is net positive and passive is net negative. Amazing finding.
6. Then subset out each of active passive and run a meta. Found that both for active and passive, vegetation is a path forward for restoration (list techniques in paper) particularly active.
7. Adding water does not generate significant effects, soil remediation is powerful active intervention and ignoring passive changes in soil is a dramatic and negative impediment to restoration.
   
8. Habitat is challenging to actively restore but plants, animals, and soils are viable and significant positive outcomes that can be restored.
   
9. Passive restoration confirms that habitat does not generally respond, plants can recover passively, but soil does not.
10. Additional ideas - test ONE key moderator such as aridity index or length of experiment?